在R中找到拟合线和数据点之间的最小差异
Finding minimum difference between fit line and data points in R
目前，我正在尝试找到一种在 R 中编写函数的方法，该函数将从应力应变曲线中找到 material 的屈服强度。我有 24 个经过拉伸测试的试样的应力应变图。
通常，屈服强度是通过取应力应变曲线的线性部分并将其偏移 0.2% 得到的。该线与原始应力应变曲线相交的位置称为 material.
的屈服强度
我可以找到图形线性部分的斜率。我遇到的问题是偏移该线并找到它与原始曲线相交的位置。
参考下图：
我的应力应变图是一组离散的数据点，因此我将通过删除一些点来将线性曲线拟合到图表的第一部分。在我有一个线性方程之后，我将抵消它的 0.2%。使用偏移方程，我会将其应用于相应的应力应变曲线。
我会求最小值的绝对值，这样我就可以得到一个很好的交点近似值。如果我不使用绝对值，那么我认为 R 会发现拟合点和数据点之间存在很大的负差异，远离线性方程离开页面的地方。
为了快速运行我的代码，请从Dropbox下载csv文件link。
#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 是我要估计的点。
谢谢！
使用末尾注释中的 stress 和 strain 计算应力-应变曲线 (spl) 的样条近似值，并使用 uniroot 计算其路口，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)
查看以下代码：
# 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])
目前，我正在尝试找到一种在 R 中编写函数的方法，该函数将从应力应变曲线中找到 material 的屈服强度。我有 24 个经过拉伸测试的试样的应力应变图。
通常，屈服强度是通过取应力应变曲线的线性部分并将其偏移 0.2% 得到的。该线与原始应力应变曲线相交的位置称为 material.
的屈服强度我可以找到图形线性部分的斜率。我遇到的问题是偏移该线并找到它与原始曲线相交的位置。
参考下图：
我的应力应变图是一组离散的数据点，因此我将通过删除一些点来将线性曲线拟合到图表的第一部分。在我有一个线性方程之后，我将抵消它的 0.2%。使用偏移方程，我会将其应用于相应的应力应变曲线。
我会求最小值的绝对值，这样我就可以得到一个很好的交点近似值。如果我不使用绝对值，那么我认为 R 会发现拟合点和数据点之间存在很大的负差异，远离线性方程离开页面的地方。
为了快速运行我的代码，请从Dropbox下载csv文件link。
#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 是我要估计的点。
谢谢！
使用末尾注释中的 stress 和 strain 计算应力-应变曲线 (spl) 的样条近似值，并使用 uniroot 计算其路口，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)
查看以下代码：
# 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])