Hyperspectral images can be seen as a generalisation

Assignment Programming for Data Science (301113) Due: 9 Oct 2017 AEDT 23:59 Centre for Research in Mathematics School of Computing, Engineering and Mathematics Question 1. Manipulate spectral image (10 points) Hyperspectral images can be seen as a generalisation of normal colour images such as RGB images. In a normal RGB colour image, there are 3 channels, i.e. channels for red colour, green colour and blue colour. It is normally organised in a 3D array, e.g. an array X of size 100×120×3, for a raster image with 100 pixels in its vertical line and 120 pixels in horizontal line. Then X[ , , 1], the first “slice” of the array which is a matrix, is the matrix containing all red colour information typically ranging from 0 to 1 with 0 meaning totally dark and 1 brightest in terms of red colour presence. Similarly X[ , , 2] is for green and X[ , , 3] for blue. Hyperspectral image (HSI) has similar structure as RGB images. Besides RGB channels, a typical HSI has hundreds other channels (sometimes called bands or wavelengths). R can render an RGB image, so in the same way, HSI can also be displayed as an image but using its RGB channels. Another way to look at HSI is that each pixel in HSI is a spectrum. For example, if X of size 610×340×103 is a 3D array storing an HSI, then there are 103 channels in this HSI and X[2,3,] is a vector of length 103 which is the spectrum at the second row and the third column (imagine an HSI as a cube). Due to the rich information contained in the bands, HSI is usually used to characterise the objects in the image. Each pixel is either a pure material or a mixture of several materials. Here we will look at this HSI called Pavia University scene remotely sensed by ROSIS sensor http://www.opairs.aero/rosis_en.html. It has 103 channels. The data can be accessed from http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes where some more information about the data set can be found. Each pixel has been classified manually as one of the 10 classes in the scene which includes a background class. The whole data set and ground truth are packed in a matlab file called paviauni.mat which is in vuws. This is the file you need to work on. You can use readMat() in R package called R.matlab to read in matlab file library(R.matlab) A <- readMat(‘paviauni.mat’) names(A) [1] “X” “groundtruth” groundtruth is the classification information of each pixel and X are the actual HSI. Write R scripts for the following tasks. 1. Produce a plot of 100 spectra randomly chosen from each class (if a class has less than 100 spectra then take them all) and plot classes in separate subfigures as the following (5 points) 0 20 40 60 80 100 0 2000 4000 0 0 20 40 60 80 100 500 1500 2500 1 0 20 40 60 80 100 0 1000 3000 2 0 20 40 60 80 100 500 1000 2000 3 0 20 40 60 80 100 0 2000 4000 4 0 20 40 60 80 100 1000 3000 5000 5 0 20 40 60 80 100 0 1000 2000 3000 6 0 20 40 60 80 100 400 800 1200 1600 7 0 20 40 60 80 100 500 1000 2000 8 0 20 40 60 80 100 0 400 800 1200 9 2. Display Pavia University data set as a colour image using channels 57, 27 and 17 as RGB channels. The function to use is grid.raster() in the package called grid. Note that the values in the image has to be normalised to the range from 0 to 1. (5 points) Question 2. Random walks on graph (10 points) Simulate random walks on a directed graph. See https://en.wikipedia.org/wiki/Directed_graph for details about directed graph. Assume there areN vertices and the connections between vertices are shown in Fig. 1. A random walker can land on any vertex (site). Then he begins to move to the next site along an arrow with the probability written on the arrow. This process repeats for e.g. 1000 times. Note that the walker can stay on the site at next step and the summuation of all outgoing probabilities of any site (arrows pointing outwards to other sites, including arrow pointing to itself ) has to be 1. See V1 in Fig. 1 as an example. Figure 1: Directed graph. Labelled circles are vertices (sites). Arrows (called edges or links) are showing possible next move. The numbers on arrows are the probabilities that the random walker can move from one site to another. Create a function for this task and a script file to demonstrate how to use it. There is no restriction on how it should be implemented. The only requirement is the final output should be a vector of the possibilities corresponding to the chance of each site to be visited based on the configuration of the directed graph. So you can design a directed graph and specify the probabilities of the move. Question 3. Cross validation for ridge regression (20 points) Cross validation (CV) is a useful tool for model selection, particularly for the selection of model hyperparameters. Some details about CV can be found in https://en.wikipedia.org/wiki/Cross-validation_(statistics). To build a model, one needs training data and test data. On training data, some models can be built. Usually, it is called model training or simply training. While the test data, is the “unseen data” not available for training, that the trained models are applied to and on which models performance indicators are calculated such as mean squared error (MSE), so that the best model can be selected. The basic idea of CV is to shuffle the data randomly first and then divide it into k partitions. Then reserve one part as test set, and take the rest as training set on which different models can be built. Each of the k parts serves as test set once so that this process is repeated k times. Therefore for one type of model, it can be trained k times, evaluated k times (on test partition) and their performance statistics can be summarised. Finally the best model can be chosen according to the statistics such as the smallest MSE. This is called the k-fold cross validation. The k-fold cross validation process can repeatN times, then it is called k-fold cross validation with repetitions. We are considering a linear regression model called ridge regression. It is similar to the simple linear regression that we have learnt before. However, ridge regression imposes a constraint on the regression coefficients to make them not “too large”. This is controlled by a hyperparameter called lambda in the model. The larger the lambda, the smaller the regression coefficients. We use linearRidge function in the ridge package for ridge regression. Download ridge package if you do not have it yet and check the help page for linearRidge function. A simple use case is shown below. library(ridge) library(MASS) head(Boston,4) # Show some observations in Boston Housing data set crim zn indus chas nox rm age dis rad tax ptratio black 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 lstat medv 1 4.98 24.0 2 9.14 21.6 3 4.03 34.7 4 2.94 33.4 ridgeregression.model <- linearRidge(medv ~ ., data=Boston, lambda=1) ridgeregression.model Call: linearRidge(formula = medv ~ ., data = Boston, lambda = 1) (Intercept) crim zn indus chas 21.023352544 -0.059891185 0.017709378 -0.072402885 2.310651531 nox rm age dis rad -3.922337411 2.875263795 -0.009292774 -0.249729427 -0.004395417 tax ptratio black lstat -0.002731648 -0.535516506 0.006194224 -0.261367653 In the above example, we use Boston data set, which is contained in MASS package. See https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/Boston.html for details. The response variable is medv component, i.e. the median value of owner-occupied homes in $1000s. Note that the hyperparameter lambda is chosen to be 1, which is totally arbitary. However, what value of lambda is reasonable? CV is the solution. You need to use k-fold CV to determine the best model in terms of MSE. MSE is defined as MSE = ΣN i=1(yi − ˆyi)2 N where yi is the ith reponse variable in test set, and ˆyi is the ith estimation of yi by using the model, which can be computed by using p
redict(ridgeregression.model,test.data) function, where ridgeregression.model is the model you trained, and test.data is the new data you reserved for testing in CV. The smaller the MSE, the better the model. Create a function to do k-fold cross validation for ridge regression using the following template: CV.RidgeRegression <- function(y,X,K=10,repetitions=1,lambdas) { # y: the response variable (should be a vector) # X: the predictors (should be a matrix) # K: the value of k in k-fold CV # repetitions: the number of repetitions in k-fold CV. # lambdas: a vector of values of lambda to be tested. # …: Any other input arguments appropriate # Your code goes from here # Return output } In CV.RidgeRegression() function, you need to provide the following functionalities: 1. Save the input data, i.e. y and X into file; 2. Save models (hyperparameters) and indices of training and test data set into file; 3. Display validation process, e.g. which fold it is running, MSE, lambda value, summary statistics of the MSE values for one lambda value (mean will be sufficient) and so on, similar to the following; MSE lambda Run 1 Fold 1 40.76059 1 Run 1 Fold 2 29.01013 1 Run 1 Fold 3 32.58421 1 Run 1 Fold 4 24.39025 1 Run 2 Fold 1 26.14238 1 Run 2 Fold 2 23.01231 1 Run 2 Fold 3 30.73557 1 Run 2 Fold 4 48.18946 1 Mean of MSE is 31.85311 for lambda = 1 MSE lambda Run 1 Fold 1 35.86531 2 Run 1 Fold 2 55.07175 2 Run 1 Fold 3 26.77029 2 Run 1 Fold 4 37.20461 2 Run 2 Fold 1 47.18468 2 Run 2 Fold 2 29.57044 2 Run 2 Fold 3 29.12745 2 Run 2 Fold 4 46.97892 2 Mean of MSE is 38.47168 for lambda = 2 The best model is when lambda = 1 4. Plot the prediction of the best model against the response with all the data with proper labels such as title, axes labels, legends and so on; 10 20 30 40 50 5 10 15 20 25 30 35 Ridge Regression with lambda = 1 Observed Median Value Prediction Model vs obs Ideal model output The code to produce the above plot is yhat <- predict(bestmodel,Boston[,1:13]) plot(Boston$medv,yhat,xlab=’Observed Median Value’,ylab=’Prediction’, main=paste(‘Ridge Regression with lambda =’,fl)) abline(coef=c(1,1)) legend(‘topleft’,c(‘Model vs obs’,’Ideal model output’),lty=c(NA,1),pch=c(1,NA)) You need to plug in the best model you selected from CV as bestmodel in the above code. 5. Return a list of things such as models generated and their MSE’s, best model which corresponds to the minimum MSE, MSE of applying the best model to the whole data set, and anything else that is relevant. Create a simple document to report on this function using RMarkdown. In this report, you need cover the following topics: 1. Explain the purpose of this function; 2. Explain how it is implemented; 3. Explain the input/output arguments and the usage of the function; 4. Explain behaviour of this function such as save data to file, screen output and so on; 5. Include a show case using Boston data set with actual output of your function and plots; Something like this cvresults <- CV.RidgeRegression(Boston$medv,Boston[setdiff(names(Boston),’medv’)], K=10,repetitions=2,lambdas=seq(0,2,1)) with all its outputs. 6. Any issues with the function. The idea for this report is to provide a guidance to users like a user manual.

Leave a Reply

Your email address will not be published. Required fields are marked *