在R中找到擬合直線和數據點之間的最小差異

Question

目前，我試圖找到一種方法來在R中編寫一個函數，該函數可以從應力應變曲線中找到材料的屈服強度。我有24個拉伸測試的應力應變圖。

通常，屈服強度是通過取應力應變曲線的線性部分並將其抵消0.2％而得到的。那條線與原始應力應變曲線相交的地方被稱爲材料的屈服強度。

我可以找到圖的線性部分的斜率。我遇到的問題是抵消該線並找到它與原始曲線相交的位置。

請參考下面的圖：

我的應力應變圖是一組離散的數據點，所以我將通過除去一些點擬合的曲線圖的第一部分的線性曲線。在我有一個線性方程後，我將抵消它的0.2％。使用偏移方程，我會將其應用於相應的應力應變曲線。

我會做的最小值的絕對值，以便我得到一個很好的相交點的近似值。如果我不使用絕對值，那麼我認爲R會發現在線性方程離開頁面的地方，擬合和數據點之間存在很大的負差。

爲了快速運行我的代碼，請從Dropbox鏈接下載csv文件。

CSV File

#Set the working directory to where you saved the CSV file and R script. 
setwd("C:/your_working_directory") 

#Read in the CSV 
test_file <- read.csv("C:/your_working_directory/test.csv", 
         header = TRUE, 
         quote="\"", 
         stringsAsFactors= TRUE, 
         strip.white = TRUE) 

#Assigns the values from the stress column to a vector 
stress <- test_file$stress[2:176] 

#Assigns the values from the strain column to a vector 
strain <- test_file$strain[2:176] 

#Plotting the stress and strain, I only inlcluded the first 175 points 
#so that you can see where the curve starts to 
#bend as shown in the example picture. 
plot(strain, stress, main='Stress v Strain', xlab='Strain (in/in)', ylab='Stress (PSI)') 

#-------------Get Slope Function Section---------------------# 

#get.slope function returns the slope of the passed in values 
get.slope<-function(stress,strain){ 
    LinearFit<- lm(stress~strain) 
    Slope <- summary(LinearFit)$coefficients[2, 1] 
} 

#calls function to fit first degree polynomial equation, 
#notice that only the first 100 points are used where the curve is 
#still fairly linear: 
modulus<-get.slope(stress[1:100],strain[1:100]) 

LinearFit <- lm(stress~strain) 
print(summary(LinearFit)) 

這是什麼劇情應該輸出：

而且其中紅色的X是，我想估計點。

謝謝！

Answer 1

使用stress，並從注在結束計算樣條逼近到應力 - 應變曲線（spl），並使用uniroot計算其相交strain，B.

# calculate B, the intersection point 

modulus <- coef(lm(stress ~ strain, subset = 1:100))[2] # slope 

spl <- splinefun(strain, stress) # spl is a function giving stress as a function of strain 
fun <- function(strain) modulus * (strain - 0.002) - spl(strain) # B[1] is root of fun 
out <- uniroot(fun, range(strain)) 
B <- c(out$root, spl(out$root)) # x and y coords of B, the intersection 

# create plot 

plot(stress ~ strain, type = "l", col = "red", 
    xaxs = "i", yaxs = "i", # axes start at 0 
    xaxt = "n", yaxt = "n", # no ticks or tick labels 
    xlab = epsilon ~ "(strain, in/in)", ylab = sigma ~ "(stress)") 

points(B[1], B[2], pch = 20, cex = 2, col = "lightblue") # blue dot at B 

segments(0.002, 0, B[1], B[2]) # AB 
segments(0, B[2], B[1], B[2], lty = 2) # horizontal dashed line segement 
segments(B[1], 0, B[1], B[2], lty = 2) # vertical dashed line segment 

# axis tick labels 
axis(1, 0.002) # mark 0.002 on X axis 
axis(1, B[1], ~ epsilon[y]) 
axis(2, B[2], ~ sigma[y], las = 1) 

# mark O, A and B 
text(B[1], B[2], "B", adj = c(-0.5, 1)) 
text(0.002, 0, "A", adj = c(1, -0.5)) 
text(0, 0, "O", adj = c(-1, -0.5)) 

注意：使用的輸入是（爲了將來的參考，這是如何可重複地提供輸入以保持所有內容獨立且容易複製到R中）：

stress <- c(113.3385421, 462.297649, 754.3743987, 873.7138964, 917.3587659, 
957.76731, 947.5992303, 962.0677743, 960.7792493, 955.4918091, 
969.4260236, 971.3544525, 965.0086849, 968.7318796, 969.209709, 
969.6165097, 969.0247115, 964.7810702, 977.008659, 975.3817792, 
974.3037574, 980.0322212, 966.6442819, 971.5307328, 975.0376129, 
984.5745059, 984.48346, 986.4060775, 1004.222656, 1003.195645, 
1040.114098, 1095.025409, 1315.958855, 1592.423499, 1872.966804, 
2152.901522, 2442.51519, 2718.843266, 3003.570806, 3271.090355, 
3549.170579, 3822.213577, 4101.431228, 4380.060633, 4648.839972, 
4922.631032, 5180.630799, 5440.223224, 5708.234808, 5951.168745, 
6197.33933, 6443.517986, 6679.028457, 6917.311952, 7159.027746, 
7387.715265, 7623.21379, 7844.363228, 8087.072456, 8318.090361, 
8537.992923, 8768.527188, 8981.64425, 9209.520427, 9422.6094, 
9624.035463, 9864.527304, 10051.87128, 10284.48604, 10497.43005, 
10717.77622, 10949.09568, 11149.25129, 11382.17503, 11584.37111, 
11795.36732, 12022.06442, 12224.53944, 12474.11263, 12677.37575, 
12886.3005, 13131.58969, 13329.81914, 13564.19282, 13787.43212, 
14004.68466, 14235.62444, 14437.42308, 14677.19785, 14903.37871, 
15135.13918, 15354.2317, 15576.63673, 15803.95985, 16018.27633, 
16248.07093, 16474.49392, 16688.64251, 16925.30236, 17137.32041, 
17358.54895, 17579.98605, 17790.69395, 18011.11179, 18228.41309, 
18436.13068, 18649.90766, 18837.8951, 19040.28747, 19216.29398, 
19414.22349, 19607.13483, 19791.43297, 19971.1885, 20159.73545, 
20348.92618, 20529.35777, 20709.53786, 20894.83492, 21071.76642, 
21243.59838, 21417.62835, 21590.91091, 21766.39697, 21929.31289, 
22094.33385, 22251.90131, 22412.12427, 22563.67302, 22700.23859, 
22844.01099, 22980.1562, 23105.12601, 23224.57722, 23338.65058, 
23449.11698, 23547.68608, 23646.64228, 23741.96956, 23831.31074, 
23910.64462, 23979.96378, 24062.2828, 24127.02416, 24186.16297, 
24259.75418, 24313.22296, 24360.72727, 24425.38727, 24463.31174, 
24517.21606, 24564.4227, 24607.92519, 24663.63008, 24689.33859, 
24743.03725, 24779.71259, 24824.05871, 24886.27435, 24890.73302, 
24951.50516, 24979.0333, 25015.69217, 25059.98954, 25077.32958, 
25141.987, 25154.45447, 25190.59419, 25236.6169, 25249.49796, 
25302.25032, 25322.93226, 25352.40432, 25388.76875) 

strain <- c(0, 4e-05, 8.5e-05, 0.00011, 0.00011, 0.000115, 1e-04, 8.5e-05, 
7.5e-05, 5e-05, 4.5e-05, 3.5e-05, 3e-05, 3e-05, 2.5e-05, 3.5e-05, 
2.5e-05, 1.5e-05, 2e-05, 1e-05, 2.5e-05, 2e-05, 2.5e-05, 2e-05, 
1.5e-05, 2.5e-05, 2e-05, 2.5e-05, 3e-05, 2.5e-05, 3.5e-05, 3.5e-05, 
6.5e-05, 9.5e-05, 0.000125, 0.00015, 0.00018, 0.00021, 0.00024, 
0.000275, 3e-04, 0.00033, 0.00036, 0.00039, 0.000425, 0.00045, 
0.000475, 0.000505, 0.00053, 0.000555, 0.00058, 6e-04, 0.00062, 
0.000645, 0.000665, 0.000685, 7e-04, 0.00072, 0.000735, 0.000755, 
0.000775, 0.000795, 0.000815, 0.000825, 0.000845, 0.00086, 0.000875, 
0.000895, 0.000915, 0.00093, 0.00095, 0.000965, 0.000985, 0.001, 
0.00102, 0.00104, 0.00106, 0.00108, 0.0011, 0.00112, 0.00114, 
0.00116, 0.001185, 0.0012, 0.001225, 0.001245, 0.00127, 0.00129, 
0.00131, 0.001335, 0.001355, 0.001385, 0.001405, 0.00143, 0.00145, 
0.001475, 0.0015, 0.001525, 0.001545, 0.00157, 0.001595, 0.001615, 
0.001645, 0.001665, 0.00169, 0.001715, 0.00174, 0.001765, 0.001785, 
0.001815, 0.001835, 0.00186, 0.001885, 0.001905, 0.001935, 0.001955, 
0.00198, 0.00201, 0.00204, 0.00207, 0.002095, 0.002125, 0.00216, 
0.00219, 0.00222, 0.002255, 0.00229, 0.00233, 0.002365, 0.00241, 
0.00245, 0.0025, 0.002545, 0.002595, 0.002645, 0.002695, 0.00275, 
0.002805, 0.00286, 0.002925, 0.00299, 0.003055, 0.003125, 0.0032, 
0.00328, 0.003365, 0.00345, 0.003535, 0.00363, 0.00373, 0.003825, 
0.003915, 0.00401, 0.004095, 0.004195, 0.004285, 0.004375, 0.00447, 
0.00456, 0.00466, 0.004755, 0.00485, 0.00495, 0.00505, 0.005145, 
0.005245, 0.005345, 0.005445, 0.005555, 0.005665, 0.005775, 0.00589, 
0.006005, 0.006115) 

Answer 2

看看下面的代碼：

# Create data for example 
strain <- seq(0, 2, by = 0.1) 
stress <- sin(strain) 

plot(strain, stress, type = "l") 

# Fit a linear model and plot the fitted values 
fit <-lm(strain ~ stress) 
lines(strain, fitted(fit)) 

# Find the distance on x-axis 
dist <- unname(-coef(fit)[1]/coef(fit)[2]) 
dist 
#[1] 0.1197057 
text(y = 0, x = dist, labels = round(dist, 2)) 


# Find point of intersection of curves 
indx <- which(diff(stress > fitted(fit))!=0) 
strain[indx] 
#[1] 0.1986693 

text(y = strain[indx], x = stress[indx], labels = strain[indx])