Skip to content

How to estimate the parameter for many variable?

9 messages · Jim Lemon, Rui Barradas, SITI AISYAH ZAKARIA

SITI AISYAH ZAKARIA
# 
Dear all,

Can I ask something about programming in marginal distribution for spatial
extreme?
I really stuck on my coding to obtain the parameter estimation for
univariate or marginal distribution for new model in spatial extreme.

I want to run my data in order to get the parameter estimation value for 25
stations in one table. But I really didn't get the idea of the correct
coding. Here I attached my coding

x <- data.matrix(Ozone_weekly2)
x
head(gev.fit)[1:4]
ti = matrix(ncol = 3, nrow = 888)
ti[,1] = seq(1, 888, 1)
ti[,2]=sin(2*pi*(ti[,1])/52)
ti[,3]=cos(2*pi*(ti[,1])/52)
for(i in 1:nrow(x))
  + { for(j in 1:ncol(x))
    + {x[i,j] = 1}}

My problem is highlighted in red color.
And if are not hesitate to all. Can someone share with me the procedure,
how can I map my data using spatial extreme.
For example:
After I finish my marginal distribution, what the next procedure. It is I
need to get the spatial independent value.

That's all
Thank you.
1 day later
Jim Lemon
# 
Hi Siti,
I think we need a bit more information to respond helpfully. I have no
idea what "Ozone_weekly2" is and Google is also ignorant. "gev.fit" is
also unknown. The name suggests that it is the output of some
regression or similar. What function produced it, and from what
library? "ti" is known as you have defined it. However, I don't know
what you want to do with it. Finally, as this is a text mailing list,
we don't get any highlighting, so the text to which you refer cannot
be identified. I can see you have a problem, but cannot offer any help
right now.

Jim

On Thu, Jul 8, 2021 at 12:06 AM SITI AISYAH ZAKARIA
<aisyahzakaria at unimap.edu.my> wrote:
Rui Barradas
# 
Hello,

Also, in the code

x <- data.matrix(Ozone_weekly)

[...omited...]

for(i in 1:nrow(x))
   + { for(j in 1:ncol(x))
     + {x[i,j] = 1}}

not only you rewrite x but the double for loop is equivalent to


x[] <- 1


courtesy R's vectorised behavior. (The square parenthesis are needed to 
keep the dimensions, the matrix form.)
And, I'm not sure but isn't

head(gev.fit)[1:4]

equivalent to

head(gev.fit, n = 4)

?

Like Jim says, we need more information, can you post Ozone_weekly2 and 
the code that produced gev.fit? But in the mean time you can revise your 
code.

Hope this helps,

Rui Barradas


?s 11:08 de 08/07/21, Jim Lemon escreveu:
SITI AISYAH ZAKARIA
# 
Dear all,

Thank you very much for the feedback.

Sorry for the lack of information about this problem.

Here, I explain again.

I use this package to run my coding.

library(ismev)
library(mgcv)
library(nlme)

The purpose of this is I want to get the value of parameter estimation
using MLE by applying the GEV distribution.

x <- data.matrix(Ozone_weekly2)                      x refers to my data
that consists of 19 variables. I will attach the data together.
x
head(gev.fit)[1:4]
ti = matrix(ncol = 3, nrow = 888)
ti[,1] = seq(1, 888, 1)
ti[,2]=sin(2*pi*(ti[,1])/52)
ti[,3]=cos(2*pi*(ti[,1])/52)



*for(i in 1:nrow(x))  + { for(j in 1:ncol(x))
the problem in here, i don't no to create the coding. i target my
output will come out in matrix that     + {x[i,j] = 1}}
                   show the parameter estimation for 19 variable which have
19 row and 3 column*


*                                                              row -- refer
to variable (station)  ; column -- refer to parameter estimation for GEV
distribution*thank you.
On Thu, 8 Jul 2021 at 18:40, Rui Barradas <ruipbarradas at sapo.pt> wrote:

            

      
      
      
    

        

      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jim Lemon
    

    

      
      #
    

  


  

      
Hi Siti,
I think some progress has been made. You have a data set with 888 rows
and 19 columns:

ow2<-read.table("Ozone_weekly2.txt",
  header=TRUE)
dim(ow2)
[1] 888  19

The values may be parts per million ozone in the atmosphere. The
columns may represent different measuring locations and my guess is
that rows are times of measurement, perhaps in weeks given the name.
If so, you might have run ismev::gev.fit() on all of the columns,
looking for estimated maxima at each location. Is this at all close to
what you are doing?

Jim

On Fri, Jul 9, 2021 at 9:59 AM SITI AISYAH ZAKARIA

        
<aisyahzakaria at unimap.edu.my> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Rui Barradas
    

    

      
      #
    

  


  

      
Hello,

The following lapply one-liner fits a GEV to each column vector, there 
is no need for the double for loop. There's also no need to create a 
data set x.


library(ismev)
library(mgcv)
library(EnvStats)

Ozone_weekly2 <- read.table("~/tmp/Ozone_weekly2.txt", header = TRUE)

# fit a GEV to each column
gev_fit_list <- lapply(Ozone_weekly2, gev.fit, show = FALSE)

# extract the parameters MLE estimates
mle_params <- t(sapply(gev_fit_list, '[[', 'mle'))

# assign column names
colnames(mle_params) <- c("location", "scale", "shape")

# see first few rows
head(mle_params)



The OP doesn't ask for plots but, here they go.


y_vals <- function(x, params){
   loc <- params[1]
   scale <- params[2]
   shape <- params[3]
   EnvStats::dgevd(x, loc, scale, shape)
}
plot_fit <- function(data, vec, verbose = FALSE){
   fit <- gev.fit(data[[vec]], show = verbose)
   x <- sort(data[[vec]])
   hist(x, freq = FALSE)
   lines(x, y_vals(x, params = fit$mle))
}

# seems a good fit
plot_fit(Ozone_weekly2, 1)       # column number
plot_fit(Ozone_weekly2, "CA01")  # col name, equivalent

# the data seems gaussian, not a good fit
plot_fit(Ozone_weekly2, 4)       # column number
plot_fit(Ozone_weekly2, "CA08")  # col name, equivalent



Hope this helps,

Rui Barradas


?s 00:59 de 09/07/21, SITI AISYAH ZAKARIA escreveu:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      SITI AISYAH ZAKARIA
    

    

      
      #
    

  


  

      
Dear Rui ang Jim,

Thank you very much.

Thank you Rui Barradas, I already tried using your coding and I'm grateful
I got the answer.

ok now, I have some condition on the location parameter which is cyclic
condition.

So, I will add another 2 variables for the column.

and the condition for location is this one.

 ti[,1] = seq(1, 888, 1)
  ti[,2]=sin(2*pi*(ti[,1])/52)
  ti[,3]=cos(2*pi*(ti[,1])/52)
  fit0<-gev.fit(x[,i], ydat = ti, mul=c(2, 3))

Thank you again.

        
On Fri, 9 Jul 2021 at 14:38, Rui Barradas <ruipbarradas at sapo.pt> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rui Barradas
    

    

      
      #
    

  


  

      
Hello,

With the condition for the location it can be estimated like the following.


fit_list2 <- gev_fit_list <- lapply(Ozone_weekly2, gev.fit, ydat = ti, 
mul = c(2, 3), show = FALSE)
mle_params2 <- t(sapply(fit_list2, '[[', 'mle'))
# assign column names
colnames(mle_params2) <- c("location", "scale", "shape", "mul2", "mul3")
head(mle_params2)


Hope this helps,

Rui Barradas

?s 09:53 de 09/07/21, SITI AISYAH ZAKARIA escreveu:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      SITI AISYAH ZAKARIA
    

    

      
      #
    

  


  

      
Dear Rui and Jim,

Thank you very much for your feedback.

Yes, now I get the output. And after this I will use this output as the
marginal distribution to continue the analysis on spatial extremes.


Thanks again. See you later.

        
On Fri, 9 Jul 2021 at 17:18, Rui Barradas <ruipbarradas at sapo.pt> wrote: