[R-sig-Geo] MLE of multivariate normal likelihood with large covariance matrix
Juan Du
dujuan at stt.msu.edu
Tue Jun 3 06:51:13 CEST 2008
Dear R helper,
I am trying to calculate the MLE of the Gaussian process with certain covariance structure, say, exponential, by using nlm. It works fine for small sample size like n<300, but when I increase it to 500 or more, the warning comes out and the estimate will always be initial value. Could you please give some hint about what I'm doing wrong or help me find another way to calculate the MLE of multivariate normal likelihood with large covariance matrix?
###########################
z<-seq(0,to=1,length=1000)
dis<-rdist(z)
V<-Exponential(dis, alpha =5, phi = 1)
llth <- function(para)
{
x<-y
dis<-dis
Vn=Exponential(dis, alpha=theta, phi=sigma2)
n=length(x)
logd=log(det(Vn))
xvx=mahalanobis(x, center = 0, cov = Vn)
(n/2)*log(2*pi)+(1/2)*logd+(1/2)*xvx
}
np=length(z)
#set.seed(2)
y=mvrnorm(mu=rep(0, np), Sigma=V) # This is the simulated data
length(y)
t1<-proc.time()
est1<-nlm(llth, p=4.5)
Warning messages:
1: In nlm(llth, p = 4.5) : NA/Inf replaced by maximum positive value
2: In nlm(llth, p = 4.5) : NA/Inf replaced by maximum positive value
##########################################
Thanks in advance,
Sincerely yours,
Juan
________________________________________
From: r-sig-geo-bounces at stat.math.ethz.ch [r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of r-sig-geo-request at stat.math.ethz.ch [r-sig-geo-request at stat.math.ethz.ch]
Sent: Thursday, May 29, 2008 6:00 AM
To: r-sig-geo at stat.math.ethz.ch
Subject: R-sig-Geo Digest, Vol 57, Issue 25
Send R-sig-Geo mailing list submissions to
r-sig-geo at stat.math.ethz.ch
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
r-sig-geo-request at stat.math.ethz.ch
You can reach the person managing the list at
r-sig-geo-owner at stat.math.ethz.ch
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."
Today's Topics:
1. Re: spatial probit (Eelke Folmer)
2. Problem overlying Kriging grid with measurement points
(Mauricio Zambrano)
3. Re: Problem overlying Kriging grid with measurement points
(Roger Bivand)
4. Re: Problem overlying Kriging grid with measurement points
(Mauricio Zambrano)
----------------------------------------------------------------------
Message: 1
Date: Wed, 28 May 2008 13:53:58 +0200
From: "Eelke Folmer" <E.O.Folmer at rug.nl>
Subject: Re: [R-sig-Geo] spatial probit
To: "'Bjarke Christensen'" <Bjarke.Christensen at sydbank.dk>,
<r-sig-geo at stat.math.ethz.ch>
Message-ID: <200805281154.m4SBsI7J015242 at hypatia.math.ethz.ch>
Content-Type: text/plain; charset="us-ascii"
Roger Bivand mentions LeSage's matlab toolbox
(http://www.spatial-econometrics.com/), including very clear manuals with
theory. I am unable to judge the theoretical validity of the routines but I
experienced that I was able to estimate parameters of spatial probit
regression models (with MCMC) that reflected the parameters used in the
data-generating process.
Regards,
Eelke Folmer
-----Original Message-----
From: Bjarke Christensen [mailto:Bjarke.Christensen at sydbank.dk]
Sent: 27 May 2008 12:41
To: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] spatial probit
Agus,
I recently asked a very similar question, and got the following thoughtful
reply from Roger Bivand:
https://stat.ethz.ch/pipermail/r-sig-geo/2008-May/003527.html
For my purpose, which is to estimate probabilities for use in an inverse
probability weighted regression, I've decided to explore alternative
routes, but if you find a way to estimate a probit model with a spatially
lagged limited dependent variable as a predictor, I'd be interested in
hearing about it as well.
Bjarke Christensen
------------------------------
Message: 2
Date: Wed, 28 May 2008 16:21:18 +0200
From: "Mauricio Zambrano" <hzambran.newsgroups at gmail.com>
Subject: [R-sig-Geo] Problem overlying Kriging grid with measurement
points
To: r-sig-geo at stat.math.ethz.ch
Message-ID:
<63d616b0805280721s3b026989g5c055d14e3ccb5c3 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Dear List,
I carried out a IDW interpolation and I got a plot of the
interpolations values over the grid that I defined.
However, for a better interpretation of the results, I would like to
plot, over the Kriging grid, the points that were used as data for the
interpolation, which are stored in a shapefile.
For doing this, I tried with:
spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
but I only get the plot of the interpolated value over the grid, and
the plot shows the following error message:
"error using packet 1"
"Missing 'data' argument, no omission value"
Could you give some hint about what I'm doing wrong ?.
Additional information is:
class(pp.idw)
[1] "SpatialPixelsDataFrame"
attr(,"package")
[1] "sp"
class(meteo_pt)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
And complete tried procedure:
p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")
meteo <- read.csv2("Meteorological_by_basin.csv")
# Selecting only those stations that have a NON "NA" value on the
field "PP_DAILY_MEAN_MM":
meteo_nna <- meteo[!is.na(meteo[,"PP_DAILY_MEAN_MM"]),]
# Settting the COORDINATES
coordinates(meteo_nna) <- ~EAST_ED50 + NORTH_ED50
# Projecting the coordinates of the meteorological stations into the
right system
proj4string(meteo_nna) = p4s
#proj4string(meteo_nna) <- CRS("+init=epsg:28992") ; # another way to
set up the coordinates
# Reading the boundary of the catchment
library(maptools)
catchment <- readShapePoly("only_catchment.shp", proj4string=p4s)
# Defining a sampling GRID of 1000mx1000m with regular spacing
# For avoiding a random grid (sampled randomly from the first cell),
and getting a fixed,
# and reproducable grid, it is neccessary to add the argument "offset
= c(0.5, 0.5)"
catchment.grid <- spsample(catchment, type="regular", cellsize=1000,
offset = c(0.5, 0.5))
# Makin possible that the grid can be used in the interpolations:
gridded(catchment.grid) <- TRUE
# Plotting the grid
spplot(catchment, sp.layout = list("sp.points", catchment.grid))
library(gstat)
# Interpolating with the INVERSE DISTANCE WEIGHTED
pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_nna, catchment.grid)
# Plotting the interpolated values
spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]")
# Reading the Meteorologicla Stations whitin the cachment
meteo_pt <- readShapePoints("meteo_by_basin.shp", proj4string=p4s)
#Attempting to get the overlay:
spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
Thanks in advance,
Mauricio
--
Linux user #454569 -- Ubuntu user #17469
------------------------------
Message: 3
Date: Wed, 28 May 2008 16:31:49 +0200 (CEST)
From: Roger Bivand <Roger.Bivand at nhh.no>
Subject: Re: [R-sig-Geo] Problem overlying Kriging grid with
measurement points
To: Mauricio Zambrano <hzambran.newsgroups at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <Pine.LNX.4.64.0805281630020.24036 at reclus.nhh.no>
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
On Wed, 28 May 2008, Mauricio Zambrano wrote:
> Dear List,
>
> I carried out a IDW interpolation and I got a plot of the
> interpolations values over the grid that I defined.
>
> However, for a better interpretation of the results, I would like to
> plot, over the Kriging grid, the points that were used as data for the
> interpolation, which are stored in a shapefile.
>
> For doing this, I tried with:
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
> list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
>
> but I only get the plot of the interpolated value over the grid, and
> the plot shows the following error message:
>
> "error using packet 1"
> "Missing 'data' argument, no omission value"
>
> Could you give some hint about what I'm doing wrong ?.
sp.layout=list("sp.points", meteo_pt)
See Note in ?spplot
Roger
>
> Additional information is:
>
> class(pp.idw)
> [1] "SpatialPixelsDataFrame"
> attr(,"package")
> [1] "sp"
>
> class(meteo_pt)
> [1] "SpatialPointsDataFrame"
> attr(,"package")
> [1] "sp"
>
>
> And complete tried procedure:
>
> p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")
>
> meteo <- read.csv2("Meteorological_by_basin.csv")
>
> # Selecting only those stations that have a NON "NA" value on the
> field "PP_DAILY_MEAN_MM":
> meteo_nna <- meteo[!is.na(meteo[,"PP_DAILY_MEAN_MM"]),]
>
> # Settting the COORDINATES
> coordinates(meteo_nna) <- ~EAST_ED50 + NORTH_ED50
>
> # Projecting the coordinates of the meteorological stations into the
> right system
> proj4string(meteo_nna) = p4s
> #proj4string(meteo_nna) <- CRS("+init=epsg:28992") ; # another way to
> set up the coordinates
>
> # Reading the boundary of the catchment
> library(maptools)
> catchment <- readShapePoly("only_catchment.shp", proj4string=p4s)
>
> # Defining a sampling GRID of 1000mx1000m with regular spacing
> # For avoiding a random grid (sampled randomly from the first cell),
> and getting a fixed,
> # and reproducable grid, it is neccessary to add the argument "offset
> = c(0.5, 0.5)"
> catchment.grid <- spsample(catchment, type="regular", cellsize=1000,
> offset = c(0.5, 0.5))
>
> # Makin possible that the grid can be used in the interpolations:
> gridded(catchment.grid) <- TRUE
>
> # Plotting the grid
> spplot(catchment, sp.layout = list("sp.points", catchment.grid))
>
> library(gstat)
> # Interpolating with the INVERSE DISTANCE WEIGHTED
> pp.idw <- krige(PP_DAILY_MEAN_MM~1, meteo_nna, catchment.grid)
>
> # Plotting the interpolated values
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]")
>
> # Reading the Meteorologicla Stations whitin the cachment
> meteo_pt <- readShapePoints("meteo_by_basin.shp", proj4string=p4s)
>
> #Attempting to get the overlay:
> spplot(pp.idw["var1.pred"], main="Daily Mean PP, [mm]", sp.layout =
> list("SpatialPointsDataFrame", meteo_pt[,"PP_DAILY_MEAN_MM"] ) )
>
>
> Thanks in advance,
>
> Mauricio
>
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
------------------------------
Message: 4
Date: Wed, 28 May 2008 16:53:12 +0200
From: "Mauricio Zambrano" <hzambran.newsgroups at gmail.com>
Subject: Re: [R-sig-Geo] Problem overlying Kriging grid with
measurement points
To: Roger.Bivand at nhh.no
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID:
<63d616b0805280753s4a2ee1c7l44407aaf7a5706c3 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Thanks a lot !. it works.
I had read the ?spplot before, but because I'm new in this topic, I
didn't figure out how to do the overlay in the right way.
My humble opinion, and thinking in newbies, I think that it could be a
good idea to add an example about this (just one line) to the next
version of the documentation of the gstat package, and/or to the
gstat_package-tutorial.pdf file.
Again, thanks a lot.
--
Mauricio
Linux user #454569 -- Ubuntu user #17469
------------------------------
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
End of R-sig-Geo Digest, Vol 57, Issue 25
More information about the R-sig-Geo
mailing list