daily pastebin goal
86%
SHARE
TWEET

Untitled

a guest Feb 23rd, 2018 62 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #Pre-class work
  2. dat=data.frame(x=c(1,2,3,4,5,6),
  3.                y=c(1,3,5,6,8,12))
  4.  
  5. min.RSS <- function(data, par) {
  6.   with(data, sum((par[1] + par[2] * x - y)^2))
  7. }
  8.  
  9. result <- optim(par = c(0, 1), min.RSS, data = dat)
  10.  
  11. library(Matching)
  12. data(lalonde)
  13.  
  14. optim(par=c(0,1),min.RSS,data=data.frame(y=lalonde$re78,x=lalonde$educ))
  15.  
  16. lm(re78~educ, data=lalonde)
  17.  
  18. 't
  19. Optimization:
  20.  
  21. $par
  22. [1] 920.8695 429.5908
  23.  
  24. $value
  25. [1] 19262169077
  26.  
  27. $counts
  28. function gradient
  29.      123       NA
  30.  
  31. $convergence
  32. [1] 0
  33.  
  34. $message
  35. NULL
  36.  
  37. Regression:
  38.  
  39. lm(formula = re78 ~ educ, data = lalonde)
  40.  
  41. Coefficients:
  42. (Intercept)         educ  
  43.       918.2        429.9
  44.  
  45. 't
  46.  
  47. claw <- function(xx) {
  48.   x <- xx[1]
  49.   y <- (0.46*(dnorm(x,-1.0,2.0/3.0) + dnorm(x,1.0,2.0/3.0)) +
  50.           (1.0/300.0)*(dnorm(x,-0.5,.01) + dnorm(x,-1.0,.01) + dnorm(x,-1.5,.01)) +
  51.           (7.0/300.0)*(dnorm(x,0.5,.07) + dnorm(x,1.0,.07) + dnorm(x,1.5,.07)))
  52.   return(y)
  53. }
  54.  
  55. #Claw
  56.  
  57. #clawx=sample((-200:200),replace=1)
  58. #clawdata=data.frame(y=sapply(clawx,claw),x=clawx)
  59. invclaw=function(x){-claw(x)}
  60. optim(par = -2, invclaw)
  61. optim(par = 0, invclaw)
  62. optim(par = 2, invclaw)
  63. #optim(par=2,claw,control=c(5,-1,0.01,1e-3,100,0,1e-8,0,10,1,5,1e7,0,10,10))
  64. optimize(claw,c(-2,2),maximum=1)
  65.  
  66. #Contaminated data
  67. set.seed(123)
  68. # Define the x values
  69. x <- c(1:12)
  70.  
  71. # Define the y values for the first 12 observations, in terms of x
  72. # y = -x + 10 + error term
  73. y <- c(-1*x[1:12] + 10 + rnorm(12, sd = 1))
  74.  
  75. lm(y~x)
  76.  
  77. min.reg <- function(par) {
  78.   loss <- mean((y - (par[1]*x + par[2]))^2)
  79.   return(loss)
  80. }
  81.  
  82. result.reg <- optim(par = c(1,1), min.reg)
  83. result.reg$par
  84. result.reg$value
  85.  
  86. # Introduce a totally different data generating process for 1 observation
  87. # (Maybe you have some small 'data contamination'): (21, 21)
  88. x <- c(x,21)
  89. y <- c(y, 21)
  90.  
  91. 't
  92. Use R to obtain the coefficients for the simple regression
  93. Use optim to reproduce those coefficents
  94. Plot the regression line
  95. lm1 <- reg(y~x)
  96. plot(x,y)
  97. abline(lm1, col = “blue”, lwd = 3)
  98. Consider what’s wrong with the regression line
  99. Use optim() to run a robust regression (robust to data contamination) that minimizes the median squared residual instead of the mean of the squared residuals. Try it with the default optim() method, and then experiment with different starting values. Then try it with the “SANN” method, repeating the process with different starting values.
  100. Using your best results, add a robust regression line to your plot using a different color than the one you used in (3) above.
  101. t'
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top