आर में बहुपद प्रतिगमन - वक्र पर अतिरिक्त बाधाओं के साथ

Question

मुझे पता है कि आर में मूल बहुपद प्रतिगमन कैसे करें। हालांकि, मैं केवल nls या lm का उपयोग कर सकता हूं ताकि अंक के साथ त्रुटि को कम किया जा सके।

यह ज्यादातर समय काम करता है, लेकिन कभी-कभी जब डेटा में मापन अंतर होता है, तो मॉडल बहुत प्रतिद्वंद्वी बन जाता है। क्या अतिरिक्त बाधाओं को जोड़ने का कोई तरीका है?

प्रतिलिपि प्रस्तुत करने योग्य उदाहरण:

मैं निम्नलिखित बना डेटा के लिए एक मॉडल (मेरा असली डेटा के समान) फिट करने के लिए करना चाहते हैं:

x <- c(0, 6, 21, 41, 49, 63, 166) 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 

सबसे पहले, यह साजिश करते हैं।

library(ggplot2) 
points <- ggplot(df, aes(x,y)) + geom_point(size=4, col='red') 
points 

ऐसा लगता है कि अगर हम एक लाइन के साथ इन बातों जुड़ा हुआ है, यह दिशा बदल जाएगा 3 बार, तो चलो इसे करने के लिए एक quartic फिटिंग की कोशिश करते हैं की तरह।

lm <- lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4)) 
quartic <- function(x) lm$coefficients[5]*x^4 + lm$coefficients[4]*x^3 + lm$coefficients[3]*x^2 + lm$coefficients[2]*x + lm$coefficients[1] 

points + stat_function(fun=quartic) 

ऐसा लगता है कि मॉडल अंक बहुत अच्छी तरह से फिट बैठता है ... सिवाय, क्योंकि हमारे डेटा 63 और 166 के बीच एक बड़ा अंतर को था, वहाँ एक बड़ी कील जो होने के लिए कोई कारण नहीं है है मॉडल में (मेरी वास्तविक डेटा के लिए मुझे पता है कि वहाँ कोई बड़ा चोटी है कि)

तो इस मामले में सवाल यह है:

• मैं कैसे निर्धारित कर सकते हैं कि स्थानीय अधिकतम (166, 9.8) पर होने की?

अगर यह संभव नहीं है, फिर एक और तरीके से करना यह होगा:

• मैं कैसे y = 9.8 से बड़ा बनने से y- मानों लाइन ने भविष्यवाणी की सीमित कर सकते हैं।

या शायद उपयोग करने के लिए एक बेहतर मॉडल है? (इसे टुकड़े के अनुसार करने के अलावा)। मेरा उद्देश्य ग्राफ के बीच मॉडल की विशेषताओं की तुलना करना है।

Answer 1

spline फ़ंक्शन का प्रकार आपके डेटा को पूरी तरह से मेल करेगा (लेकिन पूर्वानुमान उद्देश्य के लिए नहीं है)। स्पिन वक्र का व्यापक रूप से सीएडी क्षेत्रों में उपयोग किया जाता है और कभी-कभी यह गणित में डेटा बिंदु फिट बैठता है और अवसाद की तुलना में भौतिकी की कमी हो सकती है। here में और अधिक जानकारी here में एक महान पृष्ठभूमि परिचय।

example(spline) आपको बहुत सारे फैंसी उदाहरण दिखाएगा, और असल में मैं उनमें से एक का उपयोग करता हूं।

इसके अलावा, यह अधिक अधिक डेटा बिंदुओं के नमूने के लिए और उचित तब तक भविष्यवाणी के लिए lm या nls प्रतिगमन फिट हो जाएगा।

नमूना कोड:

library(splines) 

x <- c(0, 6, 21, 41, 49, 63, 166) 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 

s1 <- splinefun(x, y, method = "monoH.FC") 

plot(x, y) 
curve(s1(x), add = TRUE, col = "red", n = 1001) 

एक और दृष्टिकोण मैं सोचा बाधा को मापदंडों प्रतिगमन में की सीमा है कर सकते हैं, ताकि आप अपनी उम्मीद रेंज में भविष्यवाणी की जानकारी प्राप्त कर सकते।

नीचे optim के साथ एक बहुत ही सरल कोड, लेकिन केवल एक विकल्प है।

dat <- as.data.frame(cbind(x,y)) 
names(dat) <- c("x", "y") 

# your lm 
# lm<-lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4)) 

# define loss function, you can change to others 
min.OLS <- function(data, par) { 
     with(data, sum(( par[1]  + 
         par[2] * x + 
         par[3] * (x^2) + 
         par[4] * (x^3) + 
         par[5] * (x^4) + 
         - y)^2) 
      ) 
} 

# set upper & lower bound for your regression 
result.opt <- optim(par = c(0,0,0,0,0), 
       min.OLS, 
       data = dat, 
       lower=c(3.6,-2,-2,-2,-2), 
       upper=c(6,1,1,1,1), 
       method="L-BFGS-B" 
) 

predict.yy <- function(data, par) { 
       print(with(data, ((
        par[1]  + 
        par[2] * x + 
        par[3] * (x^2) + 
        par[4] * (x^3) + 
        par[5] * (x^4)))) 
       ) 
    } 

    plot(x, y, main="LM with constrains") 
    lines(x, predict.yy(dat, result.opt$par), col="red") 

Answer 2

मैं स्थानीय प्रतिगमन के लिए जाना होगा के रूप में eipi10 का सुझाव दिया। हालांकि, यदि आप चाहते हैं तो बहुपद रिग्रेशन होने के लिए, आप वर्गों के दंडित राशि को कम करने का प्रयास कर सकते हैं।

एक उदाहरण है जहाँ समारोह "बहुत ज्यादा" सीधी रेखा से भी घूम के लिए दंडित किया जाता है:

library(ggplot2) 
library(maxLik) 
x <- c(0, 6, 21, 41, 49, 63, 166)/100 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 
points <- ggplot(df, aes(x,y)) + geom_point(size=4, col='red') 

polyf <- function(par, x=df$x) { 
    ## the polynomial function 
    par[1]*x + par[2]*x^2 + par[3]*x^3 + par[4]*x^4 + par[5] 
} 
quarticP <- function(x) { 
    polyf(par, x) 
} 
## a evenly distributed set of points, penalize deviations on these 
grid <- seq(range(df$x)[1], range(df$x)[2], length=10) 

objectiveF <- function(par, kappa=0) { 
    ## Calculate penalized sum of squares: penalty for deviating from linear 
    ## prediction 
    PSS <- sum((df$y - polyf(par))^2) + kappa*(pred1 - polyf(par))^2 
    -PSS 
} 

## first compute linear model prediction 
res1 <- lm(y~x, data=df) 
pred1 <- predict(res1, newdata=data.frame(x=grid)) 
points <- points + geom_smooth(method='lm',formula=y~x) 
print(points) 

## non-penalized function 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0)) 
par <- coef(res) 
points <- points + stat_function(fun=quarticP, col="green") 
print(points) 

## penalty 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0), kappa=0.5) 
par <- coef(res) 
points <- points + stat_function(fun=quarticP, col="yellow") 
print(points) 

दंड 0.5 दिखता साथ परिणाम इस प्रकार है: आप दंड समायोजित कर सकते हैं, और grid, वे स्थान जहां विचलन दंडित किया जाता है।

Answer 3

ओट टूमेट्स स्रोत मेरे लिए काम नहीं करता था, कुछ त्रुटियां थीं। यहाँ एक सही संस्करण (ggplot2 का उपयोग किए बिना) है: `geom_smooth (विधि =" एल एम ", se:

library(maxLik) 
x <- c(0, 6, 21, 41, 49, 63, 166)/100 
y <- c(3.3, 4.2, 4.4, 3.6, 4.1, 6.7, 9.8) 
df <- data.frame(x, y) 

polyf <- function(par, x=df$x) { 
    ## the polynomial function 
    par[1]*x + par[2]*x^2 + par[3]*x^3 + par[4]*x^4 + par[5] 
} 
quarticP <- function(x) { 
    polyf(par, x) 
} 
## a evenly distributed set of points, penalize deviations on these 
grid <- seq(range(df$x)[1], range(df$x)[2], length=10) 

objectiveF <- function(par, kappa=0) { 
    ## Calculate penalized sum of squares: penalty for deviating from linear 
    ## prediction 
    PSS <- sum((df$y - polyf(par))^2) + kappa*(pred1 - polyf(par, x=grid))^2 
    -PSS 
} 

plot(x,y, ylim=c(0,10)) 

## first compute linear model prediction 
res1 <- lm(y~x, data=df) 
pred1 <- predict(res1, newdata=data.frame(x=grid)) 
coefs = coef(res1) 
names(coefs) = NULL 
constant = coefs[1] 
xCoefficient = coefs[2] 
par = c(xCoefficient,0,0,0,constant) 

curve(quarticP, from=0, to=2, col="black", add=T) 


## non-penalized function 
res <- maxBFGS(objectiveF, start=c(0,0,0,0,0)) 
par <- coef(res) 
curve(quarticP, from=0, to=2, col="red", add=T) 

## penalty 
res2 <- maxBFGS(objectiveF, start=c(0,0,0,0,0), kappa=0.5) 
par <- coef(res2) 
curve(quarticP, from=0, to=2, col="green", add=T)