转载 迭代重复加权最小二乘法(IRLS)_熙箔 加权最小二乘法matlab

1. Theory

http://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares(迭代重复加权最小二乘法)

http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_logistic_sect033.htm(SAS Online Help)

2. Implementation in R

The following R cod e uses IRLS algorithm to estimate parametersof Logistic Regression model, and the sample data we used is fromMachine Learning Ex4 - Logistic Regression. Itis a dataset representing 40 students who were admitted to collegeand 40 students who were not admitted, and their correspondinggrades for 2 exams. Our mission, is to build a binaryclassification model that estimates college admission chances basedon a student's scores on two exams(test1 and test2).

To run the codes in R, please follow these steps:

1. Download the csv file by clickling the URL:

http://spreadsheets.google.com/pub?key=0AnypY27pPCJydC1vRVEzM1VJQnNneFo5dWNzR1F5Umc&output=csv

2. Run the R scripts to build the model

# reading data
mydata =read.csv("c://stat/ex4.csv", header = TRUE)

# plots
plot(mydata$test1[mydata$admitted ==0], mydata$test2[mydata$admitted == 0], xlab="test1", ylab="test2",, col="red")
points(mydata$test1[mydata$admitted == 1],mydata$test2[mydata$admitted == 1], col="blue", pch=2)
legend("bottomright", c("not admitted", "admitted"), pch=c(1, 2),col=c("red", "blue") )

# build the model
g <- function(x) { return(log(x/(1 - x)))}

g.inv <- function(y) { return(exp(y)/(1 +exp(y))) }

g.prime = function(x) {
return(1/(x * (1 - x)))
}

V = function(mu) {
return(mu * (1 - mu))
}

niter = 10

mu = rep(mean(mydata$admitted),length(mydata$admitted))

for (i in 1:niter) {

Z = g(mu) + g.prime(mu) * (mydata$admitted -mu)

W = 1/(g.prime(mu)^2 * V(mu))

Z.lm = lm(Z ~ test1 + test2, weights = W, data =mydata)

eta = predict(Z.lm)

mu = g.inv(eta)

beta = Z.lm$coef

}

# print estimatedparameters

print(beta)

(Intercept) test1 test2
-16.3787434 0.14834080.1589085

# get 2 points (that will define a line)
x1 = c(min(mydata$test1), max(mydata$test1))

# when ax0 + bx2 + cx3 = 0 is the middle(decisionboundary line),
# so given x1 from sample data, solving to x2, we get:
x2 = (-1/beta[3]) * ((beta[2] * x1) + beta[1])

# plot
plot(x1,x2, type='l', xlab="test1", ylab="test2")
points(mydata$test1[mydata$admitted == 0],mydata$test2[mydata$admitted == 0], col="red")
points(mydata$test1[mydata$admitted == 1],mydata$test2[mydata$admitted == 1], col="blue", pch=2)
legend("bottomright", c("not admitted", "admitted"), pch=c(1, 2),col=c("red", "blue") )

The original R code is from

http://www.stanford.edu/class/stats191/logistic.html,and I changed it a little.

3. Implementation inSAS

proc import datafile="c:statex4.csv"
out=ex4
dbms=csv
replace;
getnames=yes;
run;

proc logistic data=ex4 des;
model admitted=test1test2/tech=fisher;
run;quit;

爱华网本文地址 » http://www.413yy.cn/a/25101014/232563.html