suppose I have variables y, x, and weights w, eg
x = randn(100) σ = 0.1 w = rand(length(x)) y = x + σ*randn(length(x))./w and I want to estimate 1. the coefficient and the intercept (something around 0, and 1) 2. and some simple estimate for the standard deviation of the error term (eg something around 0.1). What's the recommended way of doing this in Julia? GLM allows me to fit the model, but I could not figure out how to do it with weights, and I have a hard time extracting residuals etc to calculate an estimator for sigma (maybe I am doing it wrong, but residuals is not supported for the fit object). Julia 0.5-rc2, using latest DataFrames, last released GLM. To be clear, this does what I want: function myfit(y, X, w) W = Diagonal(sqrt(w)) WX = W*X Wy = W*y β = WX \ Wy ϵ = Wy-WX*β (β, sqrt(dot(ϵ, ϵ)/(size(X, 1)-size(X, 2)))) end myfit(y, hcat(ones(length(x)), x), w) but I want to learn the "standard" way of doing this. Best, Tamas -- You received this message because you are subscribed to the Google Groups "julia-stats" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
Le lundi 22 août 2016 à 13:26 +0200, Tamas Papp a écrit :
> suppose I have variables y, x, and weights w, eg > > x = randn(100) > σ = 0.1 > w = rand(length(x)) > y = x + σ*randn(length(x))./w > > and I want to estimate > > 1. the coefficient and the intercept (something around 0, and 1) > > 2. and some simple estimate for the standard deviation of the error term > (eg something around 0.1). > > What's the recommended way of doing this in Julia? GLM allows me to fit > the model, but I could not figure out how to do it with weights, and I > have a hard time extracting residuals etc to calculate an estimator for > sigma (maybe I am doing it wrong, but residuals is not supported for the > fit object). > > Julia 0.5-rc2, using latest DataFrames, last released GLM. > > To be clear, this does what I want: > > function myfit(y, X, w) > W = Diagonal(sqrt(w)) > WX = W*X > Wy = W*y > β = WX \ Wy > ϵ = Wy-WX*β > (β, sqrt(dot(ϵ, ϵ)/(size(X, 1)-size(X, 2)))) > end > > myfit(y, hcat(ones(length(x)), x), w) > > but I want to learn the "standard" way of doing this. though the code is prepared in some places to handle them. You can fit a generalized linear model with a normal distribution instead, but weights will be interpreted as case weights, while weighted least squares would correspond to inverse-variance/precision weights. This is currently underdocumented. The code to do that would be: df = DataFrame(x=x, y=y, w=w) fit(GeneralizedLinearModel, y ~ x, df, Normal(), IdentityLink(), wts=w) Contributions would be welcome of course! Regards > Best, > > Tamas > -- You received this message because you are subscribed to the Google Groups "julia-stats" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
On Mon, Aug 22 2016, Milan Bouchet-Valat wrote: > Le lundi 22 août 2016 à 13:26 +0200, Tamas Papp a écrit: >> suppose I have variables y, x, and weights w, eg >> >> x = randn(100) >> σ = 0.1 >> w = rand(length(x)) >> y = x + σ*randn(length(x))./w >> >> and I want to estimate >> >> 1. the coefficient and the intercept (something around 0, and 1) >> >> 2. and some simple estimate for the standard deviation of the error term >> (eg something around 0.1). >> >> What's the recommended way of doing this in Julia? GLM allows me to fit >> the model, but I could not figure out how to do it with weights, and I >> have a hard time extracting residuals etc to calculate an estimator for >> sigma (maybe I am doing it wrong, but residuals is not supported for the >> fit object). >> >> Julia 0.5-rc2, using latest DataFrames, last released GLM. >> >> To be clear, this does what I want: >> >> function myfit(y, X, w) >> W = Diagonal(sqrt(w)) >> WX = W*X >> Wy = W*y >> β = WX \ Wy >> ϵ = Wy-WX*β >> (β, sqrt(dot(ϵ, ϵ)/(size(X, 1)-size(X, 2)))) >> end >> >> myfit(y, hcat(ones(length(x)), x),w) >> >> but I want to learn the "standard" way of doing this. > Unfortunately, I don't think lm/LinearModel supports weights currently, > though the code is prepared in some places to handle them. > > You can fit a generalized linear model with a normal distribution > instead, but weights will be interpreted as case weights, while > weighted least squares would correspond to inverse-variance/precision > weights. This is currently underdocumented. The code to do that would > be: > > df = DataFrame(x=x, y=y, w=w) > fit(GeneralizedLinearModel, y ~ x, df, Normal(), IdentityLink(), wts=w) Thanks. A clarifying question: GLM uses the weight to multiply the likelihood for each observation, is this correct? (that's what it looks like from sources). Also, how should I go about extracting some estimate of the variance of the error term from GLM? Is that supported? Best, Tamas -- You received this message because you are subscribed to the Google Groups "julia-stats" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. |
Free forum by Nabble | Edit this page |