raman: Raman spectrum for ochre data

Description

Excerpt of data set used by Pisu et al.(2025) Innovative method for provenance studies in cultural heritage: A new algorithm based on observables from high-resolution Raman spectra of red ochre to analyze the Raman spectrum of red ochre, typically derived from iron oxides.

Usage

data("raman")

Arguments

Format

A data frame with 2048 observations on 2 variables.

x: Points along the x-axis for which the empirical density is available in y.
y: Empirical density.

Details

The data set is used in the example section to illustrate the use of function evoigt.

Examples

Run this code

data(raman)
plot(raman$x,raman$y,type="l", main="")
 
# select fourth peak
abline(v=raman$x[782],col="red")
abline(v=raman$x[1097],col="red")
datapeak = raman[(782):(1097),]
x <- datapeak[,1]
y <- datapeak[,2]
plot(x,y,type="l", main="")

# centering and scaling
offset<-451.1; A<- 27294.56; medianx <- 407.30
x= x - medianx
y = (y-offset)/A
smooth <- smooth.spline(x, y, spar = 0.5) 
ys <- predict(smooth, x)$y
 
# inverse sampling from empirical profile (x,y)
int<- 1:length(x)
breaks <-  x[- length(x)] + (x[-1] - x[-length(x)] ) /2
breaks = c( x[1]- (x[2]-x[1])/2 ,breaks, x[length(x)]+ (x[length(x)]-x[length(x)-1])/2)
samplesize=5000
bins <- sample(int, size=samplesize, prob=ys, replace=TRUE)
samplehist <- vector()
for (j in 1:length(bins)){
  samplehist[ j ] = runif(1, min= breaks[ bins[ j ] ], max= breaks[ bins[ j ] +1 ])
}
 
# estimate using Gibbs
if (FALSE) {
res <- evoigt(samplehist)
# classic parameters
# mu, sigma and gamma point estimates
classic.par<-unlist(res["posterior mean"])
classic.par
# physical parameters
# w_G = sqrt(8 log 2) *sigma; w_L = 2* gamma 
phy.par<- classic.par*c(0, sqrt(8*log(2)),  2)
names(phy.par)<-c("mu","w_G", "w_L")
phy.par
 
# plot estimated voigt versus empirical profile (original)
plot(x + medianx, y*A+offset,type="l",lwd=2,col="black",xlab="",ylab="")
lines(x+medianx,dvoigt(x,sigma=classic.par[2],gamma=classic.par[3])*A+offset,
col="red",lwd=3,lty=1)
legend("topright",legend = c("empirical profile","fitted profile"),
col = c("black", "red"),lty = c(1, 1),lwd = c(2, 3),bty = "n")
 
# R^2 
yhat <- dvoigt(x,sigma = classic.par[2],gamma = classic.par[3])*A+offset
my<-mean(ys*A+offset)
1- sum(ys*A+offset-yhat)^2/sum((ys*A+offset-my)^2)
 }