r-gregmisc-users

[R-gregmisc-users] SF.net SVN: r-gregmisc: [1058] trunk/RMCMC

[<wa...@us...>]

Revision: 1058
          http://svn.sourceforge.net/r-gregmisc/?rev=1058&view=rev
Author:   warnes
Date:     2007-02-16 15:15:42 -0800 (Fri, 16 Feb 2007)

Log Message:
-----------
Move Examples.R to tests/ and modify it to be useful.

Added Paths:
-----------
    trunk/RMCMC/tests/
    trunk/RMCMC/tests/Examples.R

Removed Paths:
-------------
    trunk/RMCMC/R/Examples.R

Deleted: trunk/RMCMC/R/Examples.R
===================================================================
--- trunk/RMCMC/R/Examples.R	2007-02-16 23:11:53 UTC (rev 1057)
+++ trunk/RMCMC/R/Examples.R	2007-02-16 23:15:42 UTC (rev 1058)
@@ -1,190 +0,0 @@
-
-#######
-# Sample abstract Model
-#######
-
-Example.Model <- model(
-    Y = 1/Z, #Z = 1/Y  fixed relationships must be expressed with
-             # parameters on the LHS and
-             # observed values on the RHS
-    Y ~ N(X%*%Beta,sigma.2),
-    Beta ~ N(0,Beta.var.2),
-    sigma.2 ~ InvGamma( V, 1)
-    )
-    # returns an object of class 'model'
-    
-X.data <- rnorm(10)
-Y.data <- rnorm(X.data*2)  # hidden variable
-Z.data <- 1/Y.data
-    
-Example.Data <- list(
-    Beta.var.2 = 90,
-    V = 30,
-    X = X.data,
-    Z = Z.data,
-    )
-
-######
-# Sample concrete Model
-######
-
-Example.Model <- model(
-    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
-             # parameters on the LHS and
-             # observed values on the RHS
-    Y ~ N(X%*%Beta,sigma.2),
-    Beta ~ N(0,Beta.var.2),
-    sigma.2 ~ InvGamma( V, 1),
-
-    data = list(
-                Beta.var.2 = 90,
-                V = 30,
-                X = X.data,
-                Z = Z.data,
-                )
-    )
-
-#########
-# Example Concrete Model with Blocking
-#########
-
-# Block on the beta
-Example.Model <- model(
-    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
-             # parameters on the LHS and
-             # observed values on the RHS
-    Y ~ N(X%*%Beta,sigma.2)
-    Beta ~ MvN(rep(0,30),diag(Beta.var.2,30))
-    sigma.2 ~ invGamma( V, 1),
-
-    data = list(
-                Beta.var.2 = 90
-                V = 30
-                X = X.data
-                Z = Z.data
-                )
-    )
-
-
-# Block on several parameters that are not part of the same distribution
-Example.Model <- model(
-    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
-             # parameters on the LHS and
-             # observed values on the RHS
-    Y ~ MvN(X%*%Beta,diag(sigma.2,30))
-    Beta[ 1:10] ~ NvN(0,diag(Beta.var.2,10))
-    Beta[11:30] ~ MvN(0,diag(Beta.var.2,20))
-    sigma.2 ~ InvGamma( V, 1),
-    
-    data = list(
-                Beta.var.2 = 90
-                V = 30
-                X = X.data
-                Z = Z.data
-                )
-    )
-
-# Block on sigma & beta
-Example.Model <- model(
-    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
-             # parameters on the LHS and
-             # observed values on the RHS
-    Y ~ N(X%*%Beta,sigma.2)
-    Beta ~ MvN(rep(0,30),diag(Beta.var.2,30))
-    sigma.2 ~ invGamma( V, 1)
-
-    <all other variables>
-    
-    data = list(
-                Beta.var.2 = 90
-                V = 30
-                X = X.data
-                Z = Z.data
-                )
-    )
-
-
-#########
-# Example proposal
-#########
-
-Proposal.Model <- model(
-    Beta <- N(Beta,30),
-    sigma.2 <- InvGamma( 25, 1)
-)
-
-# -->  proposalGenerator <- function(Beta,sigma^2) 
-#      { 
-#          Beta.new <- rnorm(Beta,30)
-#          sigma.2.mew <- rInvGamma( 25, 1)
-#          return( list(Beta=Beta.new,
-#                        sigma.2=sigma.2.new)
-#                 )
-#      }
-#
-# -->  proposalDensity <- function(Beta,sigma.2) 
-#      { 
-#          Beta.logLik <- log(dnorm(Beta,30))
-#          sigma.2.logLik <- dInvGamma( sigma.2, 25, 1)
-#          return( sum(Beta.loglik, sigma.2.loglik)
-#      }
-#
-
-# For the last example of blocking
-# -->  proposalGenerator <- list(
-#          f1 <- function(Beta, sigma.2, <all other parameters>)
-#               { 
-#                   Beta.new <- rnorm(Beta,30)
-#                   sigma.2.mew <- rInvGamma( 25, 1)
-#                   return( list(Beta=Beta.new,
-#                                sigma.2=sigma.2.new) )
-#               }
-#           f2 <- function(Beta, sigma.2, <all other parameters>)
-#               {
-#                   <...>
-#               }
-#      )
-
-########
-# Example Sampler 
-########
-
-    
-# Convert model code + data into R functions
-makeDensity.model <- function(model, data)
-{
-    # returns a function for <logDensity> = fun(current=<current state>)
-    #   <OR> a set of functions, one for each model parameter 
-    # parse model into 'tree'
-    # manipulate tree
-    # convert tree to output
-    density.function <- function(current.state) {}
-    density.function
-
-    
-}
-
-# For Gibb Sampler
-
-MakeGibbsProposal.model <- function(model)
-{
-    # returns a function for <new state> = fun(current=<current state>)
-}
-
-# For Metropolis Proposal
-
-
-# Covert proposal model to functions
-
-MakeProposalGenrator <- function(model, data)
-{
-    # returns a function for <new state> = fun(current=<current state>)
-}
-
-MakeProposalDensity <- function(model, data)
-{
-    # returns a function for <logDensity> = fun(current=<current state>, proposed=<proposaed state>)
-}
-
-
-

Added: trunk/RMCMC/tests/Examples.R
===================================================================
--- trunk/RMCMC/tests/Examples.R	                        (rev 0)
+++ trunk/RMCMC/tests/Examples.R	2007-02-16 23:15:42 UTC (rev 1058)
@@ -0,0 +1,198 @@
+library(RMCMC)
+
+#######
+# Sample abstract Model
+#######
+
+Example.Model <- model(
+    Y = 1/Z, #Z = 1/Y  fixed relationships must be expressed with
+             # parameters on the LHS and
+             # observed values on the RHS
+    Y ~ N(X%*%Beta,sigma.2),
+    Beta ~ N(0,Beta.var.2),
+    sigma.2 ~ InvGamma( V, 1)
+    )
+    # returns an object of class 'model'
+    
+X.data <- rnorm(10)
+Y.data <- rnorm(X.data*2)  # hidden variable
+Z.data <- 1/Y.data
+    
+Example.Data <- list(
+    Beta.var.2 = 90,
+    V = 30,
+    X = X.data,
+    Z = Z.data,
+    )
+
+######
+# Sample concrete Model
+######
+
+Example.Model <- model(
+    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
+             # parameters on the LHS and
+             # observed values on the RHS
+    Y ~ N(X%*%Beta,sigma.2),
+    Beta ~ N(0,Beta.var.2),
+    sigma.2 ~ InvGamma( V, 1),
+
+    data = list(
+                Beta.var.2 = 90,
+                V = 30,
+                X = X.data,
+                Z = Z.data,
+                )
+    )
+
+#########
+# Example Concrete Model with Blocking
+#########
+
+# Block on the beta
+Example.Model <- model(
+    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
+             # parameters on the LHS and
+             # observed values on the RHS
+    Y ~ N(X%*%Beta,sigma.2)
+    Beta ~ MvN(rep(0,30),diag(Beta.var.2,30))
+    sigma.2 ~ invGamma( V, 1),
+
+    data = list(
+                Beta.var.2 = 90
+                V = 30
+                X = X.data
+                Z = Z.data
+                )
+    )
+
+
+# Block on several parameters that are not part of the same distribution
+Example.Model <- model(
+    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
+             # parameters on the LHS and
+             # observed values on the RHS
+    Y ~ MvN(X%*%Beta,diag(sigma.2,30))
+    Beta[ 1:10] ~ NvN(0,diag(Beta.var.2,10))
+    Beta[11:30] ~ MvN(0,diag(Beta.var.2,20))
+    sigma.2 ~ InvGamma( V, 1),
+    
+    data = list(
+                Beta.var.2 = 90
+                V = 30
+                X = X.data
+                Z = Z.data
+                )
+    )
+
+# Block on sigma & beta
+Example.Model <- model(
+    Y = 1/Z,  #Z = 1/Y,  fixed relationships must be expressed with
+             # parameters on the LHS and
+             # observed values on the RHS
+    Y ~ N(X%*%Beta,sigma.2)
+    Beta ~ MvN(rep(0,30),diag(Beta.var.2,30))
+    sigma.2 ~ invGamma( V, 1)
+
+    <all other variables>
+    
+    data = list(
+                Beta.var.2 = 90
+                V = 30
+                X = X.data
+                Z = Z.data
+                )
+    )
+
+
+#########
+# Example proposal
+#########
+
+Proposal.Model <- model(
+                        Beta <- N(Beta,30),
+                        sigma.2 <- InvGamma( 25, 1)
+                        )
+
+# -->  proposalGenerator <- function(Beta,sigma^2) 
+#      { 
+#          Beta.new <- rnorm(Beta,30)
+#          sigma.2.mew <- rInvGamma( 25, 1)
+#          return( list(Beta=Beta.new,
+#                        sigma.2=sigma.2.new)
+#                 )
+#      }
+#
+# -->  proposalDensity <- function(Beta,sigma.2) 
+#      { 
+#          Beta.logLik <- log(dnorm(Beta,30))
+#          sigma.2.logLik <- dInvGamma( sigma.2, 25, 1)
+#          return( sum(Beta.loglik, sigma.2.loglik)
+#      }
+#
+
+# For the last example of blocking
+# -->  proposalGenerator <- list(
+#          f1 <- function(Beta, sigma.2, <all other parameters>)
+#               { 
+#                   Beta.new <- rnorm(Beta,30)
+#                   sigma.2.mew <- rInvGamma( 25, 1)
+#                   return( list(Beta=Beta.new,
+#                                sigma.2=sigma.2.new) )
+#               }
+#           f2 <- function(Beta, sigma.2, <all other parameters>)
+#               {
+#                   <...>
+#               }
+#      )
+
+########
+# Example Sampler 
+########
+
+    
+## Convert model code + data into R functions
+#makeDensity.model <- function(model, data)
+#{
+#    # returns a function for <logDensity> = fun(current=<current state>)
+#    #   <OR> a set of functions, one for each model parameter 
+#    # parse model into 'tree'
+#    # manipulate tree
+#    # convert tree to output
+#    density.function <- function(current.state) {}
+#    density.function
+
+    
+#}
+
+## For Gibb Sampler
+
+#MakeGibbsProposal.model <- function(model)
+#{
+#    # returns a function for <new state> = fun(current=<current state>)
+#}
+
+## For Metropolis Proposal
+
+
+## Covert proposal model to functions
+
+#MakeProposalGenrator <- function(model, data)
+#{
+#    # returns a function for <new state> = fun(current=<current state>)
+#}
+
+#MakeProposalDensity <- function(model, data)
+#{
+#    # returns a function for <logDensity> = fun(current=<current state>, proposed=<proposaed state>)
+#}
+
+
+
+library(MCMCpack)
+tmp <- make.density(Example.Model, Example.Data)
+tmp(Beta=rep(1,10), sigma.2=1)
+
+
+b.bb.model <- model(
+                    )


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

[R-gregmisc-users] SF.net SVN: r-gregmisc: [1080] trunk/RMCMC

[<wa...@us...>]

Revision: 1080
          http://svn.sourceforge.net/r-gregmisc/?rev=1080&view=rev
Author:   warnes
Date:     2007-03-14 10:21:41 -0700 (Wed, 14 Mar 2007)

Log Message:
-----------
Add Barrets LOH data set

Added Paths:
-----------
    trunk/RMCMC/data/
    trunk/RMCMC/data/BarrettsLOH.txt
    trunk/RMCMC/data/loh.R

Added: trunk/RMCMC/data/BarrettsLOH.txt
===================================================================
--- trunk/RMCMC/data/BarrettsLOH.txt	                        (rev 0)
+++ trunk/RMCMC/data/BarrettsLOH.txt	2007-03-14 17:21:41 UTC (rev 1080)
@@ -0,0 +1,41 @@
+	loh	n
+1p	7	17
+1q	3	15
+2p	4	17
+2q	3	18
+3p	5	15
+3q	4	15
+4p	5	15
+4q	3	19
+5p	6	16
+5q	12	15
+6p	5	18
+6q	3	19
+7p	1	18
+7q	3	19
+8p	5	19
+8q	3	21
+9p	11	17
+9q	2	16
+10p	2	12
+10q	2	17
+11p	3	18
+11q	5	18
+12p	3	19
+12q	4	19
+13q	6	14
+14q	3	12
+15q	1	16
+16p	4	19
+16q	5	16
+17p	19	19
+17q	5	21
+18p	5	15
+18q	6	13
+19p	5	20
+19q	6	16
+20p	2	17
+20q	0	8
+21p	0	7
+21q	6	18
+22q	4	15


Property changes on: trunk/RMCMC/data/BarrettsLOH.txt
___________________________________________________________________
Name: svn:executable
   + *

Added: trunk/RMCMC/data/loh.R
===================================================================
--- trunk/RMCMC/data/loh.R	                        (rev 0)
+++ trunk/RMCMC/data/loh.R	2007-03-14 17:21:41 UTC (rev 1080)
@@ -0,0 +1 @@
+loh <- read.table("BarrettsLOH.txt")


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.